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Abstract 

This paper extends the concept of generalized polarization tensors (GPTs), which 
was previously defined for inclusions with homogeneous conductivities, to inhomoge- 
neous conductivity inclusions. We begin by giving two slightly different but equivalent 
definitions of the GPTs for inhomogeneous inclusions. We then show that, as in the 
homogeneous case, the GPTs are the basic building blocks for the far-field expansion 
of the voltage in the presence of the conductivity inclusion. Relating the GPTs to the 
Neumann-to-Dirichlet (NtD) map, it follows that the full knowledge of the GPTs allows 
unique determination of the conductivity distribution. Furthermore, we show impor- 
tant properties of the the GPTs, such as symmetry and positivity, and derive bounds 
satisfied by their harmonic sums. We also compute the sensitivity of the GPTs with 
respect to changes in the conductivity distribution and propose an algorithm for recon- 
structing conductivity distributions from their GPTs. This provides a new strategy for 
solving the highly nonlinear and ill-posed inverse conductivity problem. We demon- 
strate the viability of the proposed algorithm by preforming a sensitivity analysis and 
giving some numerical examples. 
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1 Introduction 

There are several geometric and physical quantities associated with shapes such as eigen- 
values and capacities [53]. The concept of the generalized polarization tensors (GPTs) is 
one of them. The notion appears naturally when we describe the perturbation of the elec- 
trical potential due to the presence of inclusions whose material parameter (conductivity) is 
different from that of the background. 
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To mathematically introduce the concept of GPTs, we consider the conductivity problem 
in M d , d = 2, 3: 

V-(x(K d \n) + /cx(^))Vu = inR d , 

u(x) — h(x) = 0(|a:| 1_d ) as \x\ —> oo. 

Here, f2 is the inclusion embedded in M. d with a Lipschitz boundary, x(f2) (resp. x(R rf \ ft)) 
is the characteristic function of fi (resp. R d \ f2), the positive constant k is the conductivity 
of the inclusion which is supposed to be different from the background conductivity 1, h is a 
harmonic function in R d representing the background electrical potential, and the solution 
u to the problem represents the perturbed electrical potential. The perturbation u — h due 
to the presence of the conductivity inclusion 51 admits the following asymptotic expansion 
as |a;| — > oo: 

u{x)-h{x)= V t^ d <* h ( )M a p(k,n)dPT(x), (1.2) 

I"M/3|>1 

where T is the fundamental solution of the Laplacian (see, for example, [8]). The building 
blocks M a p(k,Sl) for the asymptotic expansion (|1.2[) are called the GPTs. Note that the 
GPTs M a p(k, f2) can be reconstructed from the far- field measurements of u by a least- 
squares method. A stability analysis of the reconstruction is provided in 1. On the other 
hand, it is shown in [5] that in the full-view case, the reconstruction problem of GPTs from 
boundary data has the remarkable property that low order GPTs are not affected by the 
error caused by the instability of higher-orders in the presence of measurement noise. 

The GPTs carry geometric information about the inclusion. For example, the inverse 
GPT problem holds to be true, namely, the whole set of GPTs, {M a/3 (k, Q) : \a\, \(3\ > 1}, 
determines k and f2 uniquely [5]. The leading order GPT (called the polarization tensor 
(PT)), {M a p{k, fi) : \a\, \fi\ — 1}, provides the equivalent ellipse (ellipsoid) which represents 
overall property of the inclusion [10] [19]. Moreover, there are important analytical and 
numerical studies which show that finer details of the shape can be recovered using higher- 
order GPTs [HI 0] . The GPTs even carry topology information of the inclusion [4] . It is 
also worth mentioning that an efficient algorithm for computing the GPTs is presented in 

The notion of GPTs appears in various contexts such as asymptotic models of dilute 
composites (c/. [29l [3TJ [12]), low- frequency asymptotics of waves [23], potential theory 
related to certain questions arising in hydrodynamics [33], biomedical imaging of small 
inclusions (see [5] and the references therein) , reconstructing small inclusions [551 033 1IH] , 
and shape description [3] . Recently the concept of GPTs finds another promising application 
to cloaking and electromagnetic and acoustic invisibility. It is shown that the near-cloaking 
effect of [28] can be dramatically improved by using multi-layered structures whose GPTs 
vanish up to a certain order [TTj . 

As far as we know, the GPTs have been introduced only for inclusions with homogeneous 
conductivities or layers with constant conductivities. It is the purpose of this paper to 
extend the notion of GPTs to inclusions with inhomogeneous conductivities and use this 
new concept for solving the inverse conductivity problem. We first introduce the GPTs 
for inhomogeneous inclusions and show that exactly the same kind of far-field asymptotic 
formula as (|1.2j) holds. We also prove important properties of the GPTs such as unique 
determination of Neumann-to-Dirichlet map, symmetry, and positivity. We then provide a 
sensitivity analysis of the GPTs with respect to changes in the conductivity distribution. We 
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finally propose a minimization algorithm for reconstructing an inhomogeneous conductivity 
distribution from its high-order GPTs. We carry out a resolution and stability analysis for 
this reconstruction problem in the linearized case and present numerical examples to show 
its viability. 

The paper is organized as follows. In section 2 we introduce the GPTs for inhomogeneous 
conductivity inclusions and prove that they are the building blocks of the far-field expansion 
of the potential. Section 3 is devoted to the derivation of integral representations of the 
GPTs. We also establish a relation between the GPTs and the NtD map. In section 4 
we prove important properties of symmetry and positivity of the GPTs and obtain bounds 
satisfied by their harmonic sums. In section 5 we perform a sensitivity analysis of the GPTs 
with respect to the conductivity distribution. We also show that in the linearized case, high- 
order GPTs capture high-frequency oscillations of the conductivity. In section 6, we present 
an algorithm for reconstructing inhomogeneous conductivity distributions from their high- 
order GPTs. The algorithm is based on minimizing the discrepancy between the computed 
and measured GPTs. 



2 Contracted GPTs and asymptotic expansions 

Let a be a bounded measurable function in R d , d = 2,3, such that a — 1 is compactly 
supported and 

Ai < a < A 2 (2.1) 

for positive constants Ai and A 2 . For a given harmonic function h in Mr, we consider the 
following conductivity problem: 

f V • erVw = in R d , 

I (2.2) 
[ u(x) - h(x) = 0(|x| 1 - d ) as \x\ oo. 

In this section we derive a full far-field expansion of (u — h)(x) as \x\ — > oo. In the course 
of doing so, the notion of (contracted) generalized polarization tensors (GPT) appears nat- 
urally. 

Let B be a bounded domain in M. d with a C 1 '^-boundary dB for some < r/ < 1. We 
assume that B is such that 

supp (cr - 1) C B. (2.3) 

Suppose that B contains the origin. Let H s (dB), for s e R, be the usual L 2 -Sobolev space 
and let Hfi (dB) := {<j> G H s (dB)\ j gB <j> = 0}. For s = 0, we use the notation Lq(8B). 

The Neumann-to-Dirichlet (NtD) map A a : H Q 1/2 (dB) ->■ iJ 1/2 (9B) is defined to be 



K[g] :=u\ dB , (2.4) 



where u is the solution to 



V • crVu = in B, 
^L- n ™ /)R (( v -n\ ( 2 - 5 ) 

dB 



a— = g on dB ( I u = 0) 
dv JdB 

for g G H Q 1 (dB). The operator Ai is the NtD map when a = 1. 
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Note that (|2.2I) is equivalent to 



( V • aVu = 

Au = 



in B, 

in W l \ B, 

on dB, 

u\ + = u|_ on dB, 

, — = 0(|a;| 1_d ) as — !• oo. 



du 
dv 



— a 



du 
dv 



(2.6) 



Here and throughout this paper, the subscripts ± indicate the limits from outside and inside 
B, respectively. 

Let T(x) be the fundamental solution to the Laplacian: 



T(x) 



2ir 



■In I 



47T 



d = 2, 
d = 3. 



(2.7) 



If u is the solution to (|2.2[) . then by Green's formula we have for x £ KL d \ B 



d(u-h) 



(u-h)(x)= / T(x-y) 

JdB OV 



(y)dsy 



dT(x - y) 
di 



dB OV 



dB VVy 

r \ f dT(x — y) 

(y)dsy - / — u\ + {y)dsy, 

+ JdB av v 



(u - h)\ + (y)ds y 



where the second equality holds since h is harmonic. Let g — a§^\-. Then we have 
u\qb = A CT [g] on dB. Thus we get from the transmission conditions in (|2.6p that 

(u-h)(x)= f T(x-y)g(y)ds y - [ d -^zA K[g]{y)dSy . ( 2 . 8 ) 

JdB JdB ov y 



For x eR d \B, we have 



,' dT(x--) \ 



\BB\ 



'dB 



V(x — y)ds y on dB, 



and hence 



JdB 

Thus we get from (|2~8| and (f279|) that 



^ - A g [g](y)<fa y = r(x - y)K- l K[g]{y)ds v . 



OB 



(2.9) 



(u-h)(x)= T(x-y)A^\A 1 - A a ){g](y)da v , xeM. d \B. (2.10) 

JdB 

Here we have used the fact that Ai : H 1 ^ 2 (dB) — > H^ 2 (dB) is invertible and self-adjoint: 

<Ai[0],/> ff i,2 iif -i,2 = (ff,A 1 [/]) H1/2|H _ 1/2 , V/,. g eff- 1/2 (ai?), 
with (,) H i/2 #-i/2 being the duality pair between H~ 1 / 2 {dB) and H 1 / 2 (dB). 
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Suppose that d = 2. For each positive integer n, let u c n and u s n be the solutions to 
when h(x) — r n cos n9 and h(x) — r" smnO, respectively. Let 



9n 



dv 



and g s n := a—± 
av 



dB. 



(2.11) 



Since f|2.2D is linear, it follows that if the harmonic function h admits the expansion 

oo 

h(x) = h(0) + ^r n (a c n cos n6 + a s n smnd) (2.12) 



with x = (r cos 9, r sin 6) , then we have 

du , 



g := a 



dv 1 



(a n ff„ + a n g n ), 

n=l 

and hence 

(u-h)(x) = J2 / r(x-y)(a^Ar 1 (A 1 -A ff )K](y) + <Ar 1 (A 1 -A CT )[^](y)) ( i Sy . (2.13) 



Note that r(x — y) admits the expansion 



oo 



n=l 



cos n9 x 

r 

r n y 



smnV x n . 
cos nu v H r' sin ra/. 



(2.14) 



where C is a constant, x = r x (cos 0X, sin# x ) and y = r y (cos9 y ,sm6 y ). Expansion (|2.14[) is 
valid if |x| — > oo and y G <9.B. The contracted generalized polarization tensors are defined 
as follows (see [XT]): 



mn 


= MZM 


:= / r" 1 cos A^ 1 (Ax 

JdB 


~ K)[9n](v) ds y , 


(2.15) 


M cs 

mn 


= m^M 


:= / r™cosm6 y Af 1 (A 1 

JdB 


- K)[g s n ](y) 


(2.16) 


mn 


= M% n [a] 


:= [ C sin ml?,, A^ 1 (Ax 


- Aa)[5rJ(y) 


(2.17) 


M ss 

mn 


= M% n [a] 


:= / Csinm^A^A! 

JdB 


- A °)i9n}(y) ds y . 


(2.18) 



^From (|2.13l) and (12.14[) . we get the following theorem. 

Theorem 2.1 Let u be the solution to (2. 2]) with d — 2. If h admits the expansion 112.12]) , 
then we have 

oo * oo 



27rmr" 

m— 1 n— 1 



x sin mfc 
^— ' 27rmr 1 ' 

m— 1 n— 1 



(2.19) 



which holds uniformly as \x\ — > oo. 
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In three dimensions, we can decompose harmonic functions as follows: 

oo n 

h(x) = h(0) + J2 E a mn r n Y™(6, V ) 7 (2.20) 



n— 1 m——n 



where (r, 0,ip) is the spherical coordinate of x and Y™ is the spherical harmonic function of 
degree n and of order m. Let 



av 



dB, (2.21) 



where u mn is the solution to (|2.2[) when /i(x) = r"Y^ n (0, <£>). It is well-known (see, for 
example, [35]) that 

r(a;-») = -X; E ^XT^^^W^ 7 )^. (2-22) 
e=o k=-t 

where {r,6,Lp) and (r',9',(p') are the spherical coordinates of x and y, respectively. Analo- 
gously to Theorem 1 2. 1[ the following result holds. 

Theorem 2.2 Let u be the solution to \2. 2\) with d = 3. If h admits the expansion 12.20\) . 
then we have 

oo I oo n , j 

M)W = -EEEE ra + n" 1 ^' ^ as 1x1 °°' (2 - 23) 

where the GPT M mn ^ = M mn ^[a] is defined by 

M mnU ~ [ Y i k (0',< P y n Ai 1 (A 1 -A l7 )[g mn \(r\0',<p')d<T. (2.24) 

JdB 

We emphasize that the definitions of contracted GPTs do not depend on the choice of B 
as long as (|2.3p is satisfied. This can be seen easily from (|2.19|) and (|2.23|) (see also section 



3 Integral representation of GPTs 

In this section, we provide another definition of GPTs which is based on integral equation 
formulations as in [27, 8 . Proper linear combinations of GPTs defined in this section coincide 
with the contracted GPTs defined in the previous section. 

Let N a (x,y) be the Neumann function of problem (|2.5[) . that is, for each fixed z G B, 
N a {x,y) is the solution to 

f V -aVN{;z) = -5 z {-) inB, 

\ , 1 f 

\ aVN(-, z) ■ V \ QB = -— , J N(x, z) da(x) = 0. 

Then the function u defined by 

u(x)=Af B , (T [g](x):= [ N <T (x,y)g(y)ds y , xeB (3.2) 

JOB 



(i 



is the solution to (j2.5[) , and hence 

k (T [g\{x)=N B , rT [g]{x), xedB. 
Let S B be the single layer potential on dB, namely, 



S B [(/>}(x) 



T(x-y)4>{y)ds y , x£R d . 



(3.3) 



(3.4) 



OB 



Let the boundary integral operator ICb (sometimes called the Poincare-Neumann operator) 
be defined by 

dT 



-(x - y)4>{y)ds 



'dB Ul/ y 

It is well-known that the single layer potential Sb satisfies the trace formula 

^S B [4>}\ ± = (±^I + K* B )[cf>} on as, 



(3.5) 



where K* B is the L 2 -adjoint of Kb- We recall that \1—K* B is invertible on L\{dB) if |A| > 1/2 
(see, for example, [25 l l38 l l8]). 

Identity ()2.10j) suggests that the solution u to (|2.2j) may be represented as 



u(x) 



h(x)+S B [4>](x), x£R d \B, 
Af Bi(T [iP}(x) + C, xeB 

for some densities 4> and on dB, where the constant C is given by 

1 



C = 



\dB\ 



OB 



(h + S B [4>}) ds. 



(3.6) 



(3.7) 



In view of the transmission conditions along dB in (|2.2p , (|3.3[) and (13.51) , the pair of densities 
{(f), tp) should satisfy 



-S B [<t>] + f 5b [0] ds + A a [ilj] =h- -— [ hds 
\oB\ J dB \oB\ J dB 



-{\l + K B )[ 



, dh 

= IT 
dv 



dB. 



(3.8) 



We now prove that the integral equation (|3.8p is uniquely solvable. For that, let 

1 



5 B [0] ~S B [<t>] - 



\dB\ 



Sb\4>\ ds. 



OB 



(3.9) 



Lemma 3.1 The operator A : H- l / 2 (dB) x H~ 1/2 (dB) iJ 1/2 (9S) x H~V 2 {dB) defined 
by 

S B A 



A:= 



<2 I + ,C *b) I 



(3.10) 



is invertible. 



As an immediate consequence of Lemma 13. II we obtain the following theorem. 

where the pair 



(3.11) 



Theorem 3.1 The solution u to A2.2\) can be represented in the form 
(<f>,ip) £ H- 1 l 2 {dB) x H Q 1/2 (dB) is the solution to 



A 



"0" 









h 



IB I IdB h d s 



\8B\ JOB 
Oh 

dZ\9B 



Proof of Lemma\3j\ We first recall the invertibility of S B : H- l / 2 {dB) ->■ H 1 / 2 (dB) in 
three dimensions (see, for instance, [38]). In two dimensions this result is not anymore true. 
However, using Theorem 2.26 of [8], one can show that in two dimensions there exists a 
unique O 6 L 2 (dB) such that 



OB 



Then we have 



A 



1 and Sb [0o] — on dB. 



-{\l + K* B )[<h\ 



and 



00 









■/ 


4>o 


Job 




00 J 



J 9B {I 1 + [Mda = J 9B 00 w + /cb ) [i]da = j 9B 00 



da = 1. 



(3.12) 

(3.13) 
(3.14) 



dimensional cases to determine uniquely ((f>,ip) € H 1 ^ 2 (dB) x H Q L/l! (dB) satisfying 



-1/2, 



,4 



"0" 




7" 






,9. 



(3.15) 



for (/,£?) e Hl /2 {dB) x H Q 1/2 (dB). In fact, if (/,c?) e H^ /2 (dB) x H-^ 2 (dB), then let 
C = T^gr J aB g and let (0, "0) be the solution to 



A 



"0" 









/ 



L5-C(iJ + /C B )[0o] 
It then follows from (f37T3|) and (|3T4l) that 

A 



- C0 O 




7" 






5. 



We now show that (|3 . 1 5[> is uniquely solvable for a given (/, g) g H^ 2 {dB) x iJ 1 ^ 2 (dB). 
Observe that it is equivalent to the existence and uniqueness of the solution in H 1 (M d ) to 
the problem 

r V • aVu = in B, 

Am = in R d \ B, 

du du - (3.16) 



\ o- 



dv 



du 

dv 



u\- - u\+ = f 
[ u(x) = OHx] 1 " 1 ) 



m 

g on dB, 
on dB, 



as \x\ — »• oo. 



The injectivity of A comes directly from the uniqueness of a solution to (|2.6p . Since u is 
harmonic in M. d \ B and u(x) — 0(|a;| 1_d ) as \x\ — > oo, there exists <\> G L\(dB) such that 



u(x) = S B [4>](x), x€WL d \B. 



If we set ip = cr|^|_, then 



where C = jr^y J as it|_. Note that 

c 1 



dB \OB\ 



HL1 



s B [<P}. 



We now have from (|3~T7f and (f3~19l) that 



Furthermore, we have 



3 = V> 



f = A CT [0] + C-S B [<t>] = k a ty\ - S B [<t>]- 



(3.17) 



(3.18) 



(3.19) 



(3.20) 



(3.21) 



Thus (4>, ip) satisfies (|3.15p and the proof is complete. □ 

We can now define the GPTs associated with a using the operator A. 

Definition 3.1 Let a be a bounded measurable function in M. d , d = 2,3, such that a — 1 is 
compactly supported and holds and let B be a smooth domain satisfying i2.3\) . For a 

multi-index a G N d with \a\ > 1, let (4> a ,ip a ) G H^ 1 / 2 (dB) x Hq 1 ^ 2 (8B) be the solution to 



A 



4>a 




X " \dB\ IdB ^ d s 






V ■ Vx a 



on dB. 



(3.22) 



For another multi-index f3 G N d , define the generalized polarization tensors associated with 
the conductivity distribution o~{x) by 



M af} = M aP (a) 



x^<p a (x) ds. 



(3.23) 



Definition 13.11 of the GPTs involves the domain B satisfying (|2.3p . However, we will 
show later that GPTs for a (in fact, their harmonic combinations) are independent of the 
choice of B satisfying (|2.3p . 

When \a\ — \f3\ — 1, we denote M := (M Q ^)| Q | = |^| =1 and call it the polarization tensor 
(matrix). Sometimes we write M = (Mij)f - =1 . 

For a given harmonic function h in K d , let (<p, ip) be the solution to (13. lip . Since 



h(x) = h(0) + 
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wc have 



£ 

M>1 



d a h(0) 



a l 



(3.24) 



By p. 61) the solution u to (|2.2I) can be written as 

d a h(0) 



l(x) = h(x) + Y ^fls B [0 a }{x), x£R d \B. 



|q|>1 



Using the Taylor expansion 



+00 



1/3 1 =0 



/8! 



which holds for all x such that \x\ — > 00 while y is bounded [5], we obtain the following 
theorem. 

Theorem 3.2 For a given harmonic function h in M. d , let u be the solution to 112.2]) . The 
following asymptotic formula holds uniformly as \x\ — > 00: 



i(x)-h{x)= Y, { —^d a h(o)M a0 d f3 r(x) 



*I.I/3|>1 



a\/3 



(3.25) 



There is yet another way to represent the solution to (|3.11[) . To explain it, let A e be the 
NtD map for the exterior problem: 



where u is the solution to 



A e [g] := u\ 9B 



Au = 
du 

= 9 



1 



\dB\ 



u, 



>>L! 



in R d \ B, 
on dB, 



dv 

u(x) — 0(\x\ 1 ~ d ) as \x\ — > 00. 
Let (0, ip) be the solution to (|3"TTj) . By ([3~2"0"j) . we have 



(3.26) 



dh 



1 



Oh 



-i + Kr B )[<f,] + — 



(3.27) 



On one hand, we obtain from the second identity in (|3.27[) that f gB = 0. On the other 
hand, the first identity in (|3.27[) says that 



4> = -H-SbM] 

dv 



dh 

dv 



on dB, 



and hence 



A 8 M = S B [fl- j±- I s B [<f\ + ^ 

\OB\ Jg B 



'dh 




dv 


dB_ 



(3.28) 
(3.29) 
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Moreover, 



and therefore, 



<p = <p + — s B [4>] 

dv 



dh 
dp 



OB 



A^ip] = Ax[4>] + S B [4>] + h\oB 



\dB 



on dB, 



iS B [<j>] + h). 



(3.30) 
(3.31) 



Combining (t3T29|) and ([3311 with 

K[ip] =S B [4>] +h\ dB 



\dB 



■(S B [4>]+h) 



in (j3"TT1) yields 



(A ff -A e )M = (Ax-A e ) 
(A 1 -A a )[^]=A 1 [0]. 



dh 



Thus we readily get 



4> = A^(Ai - A CT )(A CT - A 6 )" 1 ^! - A e ) 
^ = (A ff -A e )- 1 (A 1 -A e ) 



'dh 




dv 


OB 



'dh 




dv 


dB 



(3.32) 
(3.33) 



Note that by the uniqueness of a solution to problem (|3.16|) . it is easy to see that (A CT — A e ) : 
H~ 1/2 (dB) -> H^ /2 (dB) is invertible. 

Using p.32p gives a slightly different but equivalent definition of the GPTs. 

Lemma 3.2 For all a,0 G N d , M a p, defined by can be rewritten in the following 

form: 



'dx a 




dv 


dB 



ds. 



(3.34) 



M a0 (a) = / a^Ar^Ai - A ff )(A ff - A 6 )" 1 ^ - A e ) 
Job 

Formula Q3.34p shows how to get the GPTs from the NtD maps. 

4 Properties of GPTs 

In this section, we prove important properties for the GPTs. We emphasize that the har- 
monic sums of GPTs, not individual ones, play a key role. Let I and J be finite index 
sets. Harmonic sums of GPTs are Ylaei peJ a aibpM a p where ^2 ae j o a x a and ^2fj e jbpx^ 
are harmonic polynomials. 

The following lemma will be useful later. 

Lemma 4.1 Let I and J be finite index sets. Let h\(x) := 5Z ag j a a x a and h2(x) :— 
Ylpc j bpxP harmonic polynomials and let u\ be the solution to Ii2.ty) with h\{x) in the place 
ofh{x). Then, 

Y\Y'a a bpM a p = / (*■ - l)V«i • V/i 2 da. (4.1) 

aGlPEJ jRd 
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ui(x) :- 



Proof. Let ip = J2a£i a a^a and tf> — J2aei a a4>a- Then u\ is given by 

hx{x)+S B [4>]{x), xeR d \B, 

By (|3.22p . p.23p . and the integration by parts, we see 
a a bpM a p = I h 2 {x)(t>{x)ds x 

JOB 



OB 



ds B [4>] ds B [<t>\ 



dp + di 



OB 



OB 



OB 



ds x 

dS B [<i>] 



dv 



ds x 



h 2 (ip - J ds x - [ S B [4>]^-ds x 
\ dv J J dB dv 

, / , dhi \ f dh 2 , 

h 2 \W - -5— )ds x - / (A a [tp\ - hi) ~^—ds x 
V dv ) Job dv 

[ (h 2 tp - A ff [tp] -^-) ds x 
JdB \ dv J 



, du 
h 2 (T— 
dB \ dv 



dh 2 , 
u^r— I ds, 
av 



(a - l)V/i 2 • Vui da;, 

Jb 



which concludes the proof. □ 

Identity (|4. 1[) shows in particular that the definition of (harmonic combinations of) the 
GPTs given in the previous section is independent of the choice of B. 

4.1 Symmetry 

We now prove symmetry of GPTs. 

Lemma 4.2 Let I and J be finite index sets. For any harmonic coefficients {a a \a £ 1} 
and {bp\f3 £ J}, we have 



In particular, the first-order GPT, M, is symmetric. 



(4.2) 



Proof. The symmetry property (|4.2p can be easily deduced from the proof of Lemma 14.1 
However, we give here a slightly different proof. For doing so, let 



hi{x) := ^ a a x a , h 2 (x) 



I3£J 
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By (pH2"j) . we have 

^ ^ a a bpM aj3 



aei BeJ 



OB 



h,2(x)(f>(x)ds x 



no 



h 2 Af 1 (Ai - A CT )(A CT - A^-^Aj - A e ) 



"<%1 




dv 


OB 



<>B 



W12 
dv 



(Ai — A (T )(A - — A 6 ) _1 (Ai — A e ) 



Since 
we have 



(A CT - A e ) = (Ai - A 6 )- 1 + (Act - A e )- A (Ai - A„)(Ai ~ A e ) , 



'dhi 




dv 


dB 


- A e )~ 


-1 



aei pe.J 



dh2 
dv 



(Ai - A,,) 



dh 2 
dv 



(Ai — Act) (Act — A C ) (Ai — Act) 



Oh 



Since the NtD maps, Ai, A CT , and A e , are self-adjoint, we get (|4. 2[) which concludes the proof. 
□ 

4.2 Positivity and bounds 

Let h(x) = Y^aei "a 1 " be a harmonic function in Mr and u be the solution to (|2.2I) . As in 
the proof of Lemma 14.11 we have 



du 
dv 



dh\ , 

- u— I cfs 



^ a a apM a p = I [ha 
a ,8ei JdB V 

= / L- _-2 — -h—- (u-h i 

./as V dv dv 



ds 



= [ (a\Vu\ 2 - 2V(it - /i) ■ V/i - |V/i| 2 ) + / _|V(u-/i)| 2 

JR d \B 

= [ (ojV(u - h)| 2 + 2(cr - l)V(it - ft) • V/i + (cr - l)|V/i| 2 ) 

= / CT|V(u-ft) + (l-a- 1 )Vft| 2 + / (g ~ X) |Vfr| 2 - 

We can also check the following variational principle: 

V a a a B M a8 = min [ a\Vw + (1 - a^Whl 2 + / — - - I Vh\ 2 . (4.3) 
^ wem(R*)J Rd Jb ° 

Following the same lines of proof as in [5] for the homogeneous case, we have the following 
bounds for GPTs. 
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Theorem 4.1 Let I be a finite index set. Let {a a \a € 1} be the set of coefficients such that 
h(x) := ^2 aeI a a x a is a harmonic function. Then we have 



[ — — |V/i| 2 < V a a agM a g< [ (a-l)\Vh\ 

Jb & rzi, Jb 



(4.4) 



Proof. The bound on the left-hand side is obvious since 

min / a\Vw + (I - a^Whl 2 > 0. 

By taking w = 0, we get 

£ a a a p M aP < [ i^l!|V/ l | 2 + f ^1|V^| 2 = / (a-l)\Vh\ 2 , 

which concludes the proof. □ 

The above theorem shows that if a is strictly lager than 1 then the GPTs are positive 
definite, and they are negative definite if < a < 1. Note that optimal bounds on the 
first-order GPT have been derived in [21] [30] . 



4.3 GPTs and contracted GPTs 

The contracted GPTs appeared in the asymptotic expansions as in (|2 . 19[) and (|2 . 23[) while 
the GPTs appeared in (|3.25[) . By comparing those asymptotic formulas, we obtain the 
following lemma which relates both quantities. 

Lemma 4.3 (i) If r n cos n6 = X)| Q |=n a a a;Q an d r ™ smn ^ — ^l\ a \=n a a xa * n ^ wo dimen- 
sions, then 

M rnn = £ "a^^. M mn = £ a a a p M o/3, 
\a\—m,\0\—n \a\—m,\0\—n 

\a\=m,\0\— n \a.\—m : \0\—n 



(ii) Ifr n Y?(9,tp) = '£\a\= n < nxa 



in three dimensions, then 



M mnkl = £ a™ n a k B l M a3 . 

\a\=n,\fi\=l 

Conversely, any harmonic combination of the GPTs can be recovered from the contracted 
GPTs. 



4.4 Determination of NtD map 

It is proved in [B] (see also [5J Theorem 4.9]) that the full set of harmonic combinations of 
GPTs associated with a homogeneous inclusion determines the NtD map on the boundary of 
any domain enclosing the inclusion, and hence the inclusion. In the case of inhomogeneous 
conductivity inclusions, the same proof can be easily adapted to obtain the following result. 
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Theorem 4.2 Let I and J be finite index sets. Let ai, i = 1,2, be two conductivity distri- 
butions with supp (ai — 1) C B and satisfying A2.1\) . If 

^ ad>pM a p(ai) = X! X! a abpM al3 (a 2 ) (4.5) 
aeipej aei pej 

for any harmonic coefficients a a and bp, then 

A (Jl = A CT2 o?i <9B. (4.6) 

Using uniqueness results of the Calderon problems (for example [37J HI]) one can deduce 
from (|4.6p that a\ = a 2 under some regularity assumptions on the conductivities imposed 
in those results. In two dimensions, uniqueness holds for conductivities in L°° [14] , 

5 Sensitivity analysis for GPTs 

We now consider the sensitivity of the GPTs with respect to changes in the conductivity 
distribution. Again, we suppose that a — 1 is compactly supported in a domain B. The 
perturbation of the conductivity a is given by a + £7, where e is a small positive parameter, 
7 is compactly supported in B and refers to the direction of the changes. The aim of this 
section is to derive an asymptotic formula, as e — > 0, for the perturbation 

Am : = Y Y aab ' 3 ( M ^( a + £ 7) ~ M afj {a)) , (5.1) 

where {a a } and {bp} are harmonic coefficients and / and J are finite index sets. 
Let, as above, hi and h 2 be the harmonic functions given by 

hi{x)=^a a x a , h 2 (x) = ^bpx?. 

a p 

By (I3.32|) and a direct calculation we obtain 

Am = I h 2 Ar x ((Ai - A CT+e7 )(A ff+e7 - A 6 )- 1 - (A x - A a ){A a - A 6 )" 1 ) (A x - A e )[V/n • v\ 

JdB v ' 

= [ /i 2 A- 1 ((A CT -A CT+e7 )(A CT+e7 ~A e )- 1 

JdB K 

+ (Ai - A CT ) [(A CT+e7 - A 13 )- 1 - (A CT - A e ) _1 ] ) (Ai - A e )[V/n • v\ 
= f hi A^ 1 (A 1 - A e )(A ff - A 6 )- 1 ^ - A a+n )(A a+n - A e )-\A X - A e )[V/n • v\ 

JdB 

= / (A CT - A a+£J )[g 2 ] g e 1: 

JdB 

where 

ff J = (A <H . er -A e )- 1 (A 1 -A e )[v7»i-i/] and 32 = (A CT — A e ) _1 (Ai — A e )[\7h 2 ■ v\. (5.2) 
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Since A CT is self- adjoint, we have 

Aai = ^ / (A CT - A CT+£7 )[<? 2 + (<? 2 +g{)- \ I (A CT - A (T+e7 )[ 52 ] ,g 2 
z Job z Job 

- 5 / (A ff -A a+e7 )^]^. 

z JOB 

We need the following two lemmas. 

Lemma 5.1 If ui and u 2 are the solutions of V • (ciVui) = and V • (c 2 Vu 2 ) = with 
the Neumann boundary conditions 01^7 = g and 2 ^f = g on dB, respectively, then the 
following identity holds 

[ (A CT2 - A cri )[g]gds = - /(ctj -ct 2 )(|V(ui - u 2 )\ 2 + \\7 Ul \ 2 + \\7u 2 \ 2 )dx. (5.3) 
Job 1 J b v ' 

Proof. The following identity is well-known (see, for instance, [5]): 

/ ai\V(ui - u 2 )\ 2 dx + / (<ri - cr 2 )|Vui| 2 dx = / (A CT2 - A CTl )[g] gds. 

J B J B JOB 

We also have 

/ <r 2 |V(ui - u 2 )\ 2 dx - / (in - cr 2 )|Vu 2 | 2 dx = / (A CT1 - A CT2 )[g] g«£s. 
is Jb Job 

Subtracting those two equalities we obtain (15.3[) . □ 
Lemma 5.2 There is a constant C such that 

[|A CT1 -A CT2 1| <C||c7i- ff 2 ||£co (B) . (5.4) 

Proof. Let ui and u 2 be the solutions to V-(ctiV«i) = and V-(cr 2 Vu 2 ) = with boundary 
conditions 01^7 — g and 027^7 =3011 dB, respectively, and let h be the harmonic function 
with §^ = q on 9_B. Then we have 

[ (Ji|Vui| 2 = / V7i-V«i<^- / \\7h\ 2 dx+^ f \Vui\ 2 dx, 

./ B ^ B J B <J B 

for any e > 0. Choosing e = inf B oi := ci we get 

1 



\VuA 2 dx < 



<?i 7b 



|V/i| 2 dx. 



Similarly, we get 



where o 2 :— inf B a 2 . It then follows from (|5.3|) that 

3 



/ |Vu 2 | 2 da; < / |V/i| 2 rfx, 
0"2 Jb 



/( 



(A CT2 - h ai )[g)gds 



^ |Vtti| 2 d:E + 


y ivm 2 i 2 ^ 







+ |V/l| 2 ^) |kl-(7 2 || i( » (B) 

^ CTl £2_ JB 

< C||5 , lll/-i/2(0_B) - 02 || L°° (B) ■ 
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Thus, we obtain pT3|) . □ 

With the notation (|5 .2[) in hand, let 



9l := (A CT - A 6 )" 1 ^! - A e )[V/n • f] , (5.5) 

and let u,, for i = 1, 2, be the solution to 

f V • (ctVmj) = in B, 

du *-„ «o ( 5 - 6 ) 



cr— = 3^ on SB 

ou 



Let uf, and v\ be the solutions to 



V • (crVuf) =0 in B, 

du\ e (5-7) 
cr— = g x on dB, 

V • ((ct + ej)\7v() = in £, 

Q v e (5.8) 

a— = ffl on 8B, 



and 

C V • ((cr + e7)Vv|) =0 in B 



9u| „ (5.9) 



Then by Lemma |5. II we have 



A M = \J b e 7 (|V(u 2 + u\ -v\- v\)\ 2 + |V(u 2 + u\)\ 2 + |V(«| + wj)| 2 )da; 

~ I / £7 ( |V(M2 ~ W 2)| 2 + l V "2| 2 + |V^| 2 )dx 

-\J ct(|V(«! - ^)| 2 + iVull 2 + |V<)cix. 
Lemma 15.21 yields 

IK -«i = 0(e) (5.10) 

and 

H-«illi>(B) = 0(c), j = 1,2. (5.11) 

Thus we get 

A M = \ J 7(|V(«i + u 2 )| 2 - |V Ul | 2 - |V W2 | 2 )dx + 0(e 2 ), 
to arrive at the following theorem. 

Theorem 5.1 Let I and J be finite index sets. Let u\ and m 2 be the solutions to 15. 6\) . 
Then we have 

a a bfsM afs (a + e 7 ) = ^ ^ a a bpM af) {a) + e / 7 Vui ■ Vu 2 cfe + (9(e 2 ). (5.12) 
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6 Reconstruction of an inhomogeneous conductivity dis- 
tribution 



Over the last decades, a considerable amount of work has been dedicated to the inverse 
conductivity problem. We refer, for instance, to jTTJ [35] and the references therein. 

Here, our approach is completely different. We stably recover some important features 
of inhomogeneous conductivities using their GPTs. It should be emphasized that the GPTs 
can be obtained from boundary measurements by solving a least-squares problem [I]. The 
purpose of this section is to illustrate numerically the viability of this finding. 

For doing so, we use a least-square approach (see, for instance, [M])- Let a* be the exact 
(target) conductivity (in two dimensions) and let y mn '■— M mn (a*) (omitting for the sake 
of simplicity c and s for the superscripts in contracted GPTs). The general approach is to 
minimize over bounded conductivities a the discrepancy functional 

S ^ = \ ^2 u mn \\y mn - M mn (a)\\ 2 (6.1) 

m+n<N 

for some finite number ./V and some well-chosen weights Lo mn . The weights u> mn are used 
to enhance resolved features of the conductivity as done in [3j [18]. We solve the above 
minimization problem using the gradient descent (Landweber) method. 

6.1 Frechet derivative and an optimization procedure 

Let, again for the sake of simplicity, M mn (a) = M^ n (a) be the contracted GPTs for a given 
conductivity a. The Frechet derivative in the direction of 7, M^^cr)^], is define to be 

^From (|5.12p we obtain that 

Ml„MH= / lVu n -Vu m dx, (6.2) 

J B 

where u n and 

Urn 3X6 the solutions of 

V • (itVu) =0 in B, 

811 (6-3) 
o— = (A CT -A e )- 1 (A 1 -A e )[V/i-i/l on 35, 
ov 

with h = r n cosnO and h = r m cos m0, respectively. Note that if M mn (a) is one of the other 
contracted GPTs, then h should be changed accordingly. 

One can easily see that the adjoint M' mn (<j)* of M' rnn (<j) is given by 

lC(d)*[c] = cVu m ■ Vu n , c e M. (6.4) 

The gradient descent procedure to solve the least-square problem (|6.1[) reads 

Cfc+i = dk + ^ LUmn M' mn (t7 fe )* [y mn - M mn {a k )]. (6.5) 
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In the numerical implementation, the GPTs for the exact conductivity distribution can be 
computed by using the following formula: 

M mn (a) = [ (hi- ui)a^-ds = [ hia^-ds - [ ^-u 2 ds, (6.6) 
JdB 9f J 9B dv J dB dv 

where u\ and u 2 are the solutions to 



V • (o-(z)Vui) = in B 

dui dl 
dv dv 



9ui dhi f (6-7) 

a—— — — — on dB ( Ui = 0), 



and 

f V • a\7u 2 = 



u 2 (x)-h 2 (x) =0(\x\ l - d ) \x\^oo 



(6.8) 



respectively. Here, hi — r n cos n9 and h 2 ~ r m cosmO in two dimensions. 

On the other hand, in order to compute M' mn (a)* we need to invert the operator A a — A e . 
This can be done iteratively. In fact, the least-square solution to 

(A ff -A e )[ 5 ] = /, 

is given by 

9k+i =. 9fc +c.(A CT -A e )(/-(A ff -A e )[ fffe ]), (6.9) 

where w is a positive step-size. 

In order to stably and accurately reconstruct the conductivity distribution, we use a 
recursive approach proposed in [13] (see also [3l [4j [TSj US] ) . We first minimize the discrepancy 
between the first contracted GPTs for 1 < m,n < I. Then we use the result as an initial 
guess for the minimization between the GPTs for 1 < m, n < I + 1. This corresponds to 
choosing appropriately the weights ui mn in (|6.ip . Moreover, we refine the mesh used to 
compute the reconstructed conductivity distribution every time we increase the number of 
used contracted GPTs in the discrepancy functional. 



6.2 Resolution analysis in the linearized case 

Let d = 2 and let B be a disk centered at the origin. Consider the linearized case by 
assuming that the conductivity a is given by a = k + 57, where k =^ 1 is a positive constant 
and e is a small parameter. In that case, (|6.3p can be solved explicitly. Using Theorem 15. II 
together with Lemma [4~3l one can easily see that 

( / -i(r 1 6)r m+n ~ 2 e l{m+n)0 d6 dr ) 

\JB J l<m,n<N 

can be obtained from the contracted GPTs, M mn , for 1 < m, n < N. Therefore, the higher 
is N, the better is the angular resolution in reconstructing 7. On the other hand, it is 
clear that variations of 7 that are orthogonal (in the L 2 sense) to the set of polynomials 
(r m+n_2 )i<m,n<JV can not be reconstructed from the contracted GPTs M mn , 1 < to, n < N. 
Moreover, the reconstruction of 7 near the origin (r — 0) is more sensitive to noise than 
near the boundary of B. This is in accordance with [5] 134]. 
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6.3 Numerical illustration 



In this section, for simplicity we only consider the reconstruction from contracted GPTs of a 
conductivity distribution which is radially symmetric. Many recent works have been devoted 
to the reconstruction of radially symmetric conductivities. See, for instance, P^l 1521 155] . 
Here we consider the following conductivity distribution: 

a = (0.3r 2 + 0.5r 3 + 6(r 2 - 0.5) 2 + 3.0)/3.0, (6.10) 

and apply our original approach for recovering a from the contracted GPTs M mn , for m,n < 
N. Since the conductivity distribution a is radially symmetric we have 

M™ n = M™ n =0 for all m,n, 
C = M mn = if m ? n, 

and M m := M£f m = M^m- We use M\ to estimate the constant conductivity which has the 
same first-order GPT as follows: 

ao ■= pfmT ( } 

Then we use ao as an initial guess and apply the recursive approach described below. 

Let fc* be the last iteration step, and let Em and e a be discrepancies of GPTs and the 
conductivities, i.e., 

£ M := ^2(yn-M n (a K )) 2 , y n :=M n (<r), (6.12) 

n<N 



(N represents the number of GPTs used) and 



2 



(6.13) 



Figure [T] shows the reconstructed conductivity distribution using contracted GPTs with 
N ~ 6. In this reconstruction, the errors Em and e a are given by 

e m = 9.83318e - 005, E a = 3.5043e - 005, 

after 1598 iterations. It should be noted that the conductivity is better reconstructed near 
the boundary of the inclusion than inside the inclusion itself. Figure [2] shows how fast e m 
decreases as the iteration proceeds. The sudden jump in the figure happens when we switch 
the number of GPTs from N to N + 1. Figure [3] is for the convergence history of e CT . 



7 Conclusion 

In this paper we have introduced for the first time the notion of GPTs for inhomogeneous 
conductivity inclusions. The GPTs carry out overall properties of the conductivity distri- 
bution. They can be determined from the NtD map. We have established positivity and 
symmetry properties for the GPTs. We have also analyzed their sensitivity with respect 
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Figure 1: Reconstructed conductivity distribution. 

to small changes in the conductivity. We have proposed a recursive algorithm for recon- 
structing the conductivity from the GPTs and presented a numerical example to show that 
radially symmetric conductivities can be accurately reconstructed from the GPTs. A nu- 
merical study of the use of the GPTs for solving the inverse conductivity problem will be 
the subject of a forthcoming work. A stability and resolution analysis will be performed. It 
would also be very interesting to extend the ideas of this paper to the inverse wave medium 
problems. 
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